%===============Boyi Li 20221104================%
%=========MY==ValoMC+Model+K-Wave+2DTR==========%
%==========For multiangle solution =============%

% Create the k-Wave grid
Nx = 300;           % number of grid points in the x (row) direction
Ny = 600;           % number of grid points in the y (column) direction
dx = 0.1e-3;        % grid point spacing in the x direction [m]
dy = 0.1e-3;        % grid point spacing in the y direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy);

% Create two internal structures using makeDisk

% bone=zeros(Nx,Ny);
% bone(31:76,40:42)=1;
dim=2; % [mm]
discs = makeDisc(Nx, Ny, 150, 50, dim) +makeDisc(Nx, Ny, 100, 50, dim)+makeDisc(Nx, Ny, 200, 50, dim)+makeDisc(Nx, Ny, 150, 100, dim) +makeDisc(Nx, Ny, 150, 150, dim) +makeDisc(Nx, Ny, 150, 200, dim) + makeDisc(Nx, Ny, 150, 250, dim) + makeDisc(Nx, Ny, 150, 300, dim) +makeDisc(Nx, Ny, 100, 300, dim)+makeDisc(Nx, Ny, 200, 300, dim) + makeDisc(Nx, Ny, 150, 350, dim)+ makeDisc(Nx, Ny, 150, 400, dim);




figure,

imagesc(flip(rot90(discs)));
colormap(getColorMap);
ylabel('Axial-Dirction [mm]');
xlabel('Lateral-Dirction [mm]');
c = colorbar;  % create a colorbar
colorbar;
title('Phantom');
 

figure,

imagesc((discs));

ylabel('Lateral-Dirction[mm]');
xlabel('Axial-Dirction [mm]');
c = colorbar;  % create a colorbar
colorbar;
title('Phantom');



% Define the acoustic properties

%globe
medium.sound_speed = 1500*ones(Nx, Ny);    % [m/s]
medium.density = 1000*ones(Nx, Ny);        % [kg/m^3]
%bone
medium.sound_speed(find(bone==1)) = 2500;   % [m/s]
medium.density(find(bone==1)) = 1300;   % [kg/m^3]
%Scatter
medium.sound_speed(find(discs==1)) = 2500;   % [m/s]
medium.density(find(discs==1)) = 1300;   % [kg/m^3]



vmcmesh = createGridMesh(kgrid.y_vec*1e3, kgrid.x_vec*1e3); % [m to mm]

vmcmedium.absorption_coefficient = 0.01*ones(Nx, Ny);
vmcmedium.scattering_coefficient = 0.01*ones(Nx, Ny);
vmcmedium.scattering_anisotropy = 0.9*ones(Nx, Ny);       % scattering anisotropy parameter [unitless]
vmcmedium.refractive_index = 1.3*ones(Nx, Ny);            % refractive index [unitless]
% imagesc(discs)

% Define the acoustic properties
disc_indices = find(discs==1);

vmcmedium.absorption_coefficient(disc_indices) = 0.91;
vmcmedium.scattering_coefficient(disc_indices)= 0.1;
vmcmedium.scattering_anisotropy(disc_indices) = 0.7;           % scattering anisotropy parameter [unitless]
vmcmedium.refractive_index(disc_indices) = 1.6;    % refractive index [unitless]

% Define the Gruneisen parameter describing photoacoustic efficiency
vmcmedium.gruneisen_parameter = 0.02*ones(Nx, Ny);      % [unitless]
% Set a light source with a width of 2 mm and cosinic directional profile
% in -x direction



line1_start = [-300 -150];
line1_end = [0 0];
line1_width = 4;

line2_start = [-300 150];
line2_end = [0 0];
line2_width = 4;



lightsource1 = findBoundaries(vmcmesh, 'direction', ...
                              line1_start, ...
                              line1_end,  ...
                              line1_width);

lightsource2 = findBoundaries(vmcmesh, 'direction', ...
                              line2_start, ...
                              line2_end,  ...
                              line2_width);

% 2: A Gaussian light source
%
% Create a light source with a Gaussian directivity profile. The initial
% angles with respect to a given direction (by default, normal of the
% boundary element) follow a Gaussian with sigma = 0.1

vmcboundary.lightsource(lightsource1) = {'gaussian'};
vmcboundary.lightsource_gaussian_sigma(lightsource1) = 0.1;

% Tilt the lightsource by 22.5 degrees. This time, the direction is given in
% the coordinate system of the boundary element: (0, 1) is the normal
% direction and (1, 0) is directed along the boundary element

OMG=pi/10;%pi/4;pi/5;pi/6;pi/7;pi/8;pi/9;pi/10
OMGAngle=OMG/pi*180;   
%pi/3;%pi/4;pi/5;pi/6;pi/7;pi/8;pi/9;pi/10
%60.0;45.0;36.0;30.0;25.7;22.5;20.0;18.0

vmcboundary.lightsource_direction(lightsource1,1) = sin(OMG);
vmcboundary.lightsource_direction(lightsource1,2) = cos(OMG);

% this direction was given with respect to the surface normal
vmcboundary.lightsource_direction_type(lightsource1) = {'relative'};


% 2: A Gaussian light source
%
% Create a light source with a Gaussian directivity profile. The initial
% angles with respect to a given direction (by default, normal of the
% boundary element) follow a Gaussian with sigma = 0.1

vmcboundary.lightsource(lightsource2) = {'gaussian'};
vmcboundary.lightsource_gaussian_sigma(lightsource2) = 0.1;

% Tilt the lightsource by 22.5 degrees. This time, the direction is given in
% the coordinate system of the boundary element: (0, 1) is the normal
% direction and (1, 0) is directed along the boundary element
vmcboundary.lightsource_direction(lightsource2,1) = -sin(OMG);
vmcboundary.lightsource_direction(lightsource2,2) = cos(OMG);

% this direction was given with respect to the surface normal
vmcboundary.lightsource_direction_type(lightsource2) = {'relative'};

vmcoptions.photon_count = 10e6; 




solution = ValoMC(vmcmesh, vmcmedium, vmcboundary);
figure,
patch('Faces',vmcmesh.H,'Vertices',vmcmesh.r,'FaceVertexCData', solution.element_fluence, 'FaceColor', 'flat','EdgeColor','none' ,'LineWidth',0.5);

hold on;

ylabel('Axial-Dirction [mm]');
xlabel('Lateral-Dirction [mm]');


title(num2str(OMGAngle,'%.2f'),' Angle');
% text(0, Ny/2+2, 'gaussian tilted by deg', 'HorizontalAlignment', 'center');

c = colorbar;
c.Label.String = 'Fluence [J/mm^2]';
hold off

colormap default;



